#pragma once
extern "C" {
#include <f2c.h>
#include <clapack.h>
}
#include "LatoolException.h"
#include "IDoubleMatrix.h"

namespace LatoolNet {
	using namespace System;

	ref class DoubleGeneralTridiagonalMatrix : public IDoubleMatrix {
	private:
		double * d;
		double * dl;
		double * du;
		double * du2;
		long * ipiv;
		int m_rownum;
		int m_colnum;
		bool m_isFactorized;

		DoubleGeneralTridiagonalMatrix(DoubleGeneralTridiagonalMatrix ^ orig) {
			m_rownum = orig->m_rownum;
			m_colnum = orig->m_colnum;

			d = new double[m_rownum];
			for (int i = 0; i < m_rownum; i++) {
				d[i] = orig->d[i];
			}
			du = new double[m_rownum - 1];
			dl = new double[m_rownum - 1];
			for (int i = 0; i < m_rownum - 1; i++) {
				du[i] = orig->du[i];
				dl[i] = orig->dl[i];
			}

			du2 = new double[m_rownum -2];
			for (int i = 0; i < m_rownum - 2; i++) {
				du2[i] = orig->du2[i];
			}
			ipiv = new long[m_rownum];
			for (int i = 0; i < m_rownum; i++) {
				ipiv[i] = orig->ipiv[i];
			}

			m_isFactorized = orig->m_isFactorized;
		};
	public:
		DoubleGeneralTridiagonalMatrix(int rownum, int colnum) {
			m_rownum = rownum;
			m_colnum = colnum;

			d = new double[rownum];
			for (int i = 0; i < rownum; i++) {
				d[i] = 0.0;
			}
			du = new double[rownum - 1];
			dl = new double[rownum - 1];
			for (int i = 0; i < rownum - 1; i++) {
				du[i] = 0.0;
				dl[i] = 0.0;
			}

			du2 = new double[rownum -2];
			ipiv = new long[rownum];

			m_isFactorized = false;
		};

		~DoubleGeneralTridiagonalMatrix(){
			delete[] d;
			delete[] du;
			delete[] dl;
			delete[] du2;
			delete[] ipiv;
		};

		virtual IDoubleMatrix ^ Clone() {		
			return (IDoubleMatrix^) gcnew DoubleGeneralTridiagonalMatrix(this);
		};

		virtual property int RowNum {
			int get() {
				return m_rownum;
			}
		};
		virtual property int ColNum {
			int get() {
				return m_colnum;
			}
		};
		virtual property double default[int, int] {
			double get(int i, int j) {
				if (abs(i - j) > 1) {
					return 0.0;
				}
				if (i == j) {
					return d[i];
				} else if (j - i == 1) {
					return du[i];
				} else if (i - j == 1) {
					return dl[j];
				} else {
					return 0.0; //Unreachable. 
				}
			}
			void set(int i, int j, double value) {
				if (abs(i - j) > 1) {
					if (value == 0.0) {
						return;
					} else {
						throw gcnew LatoolException("The specified element doesn't fit to tridiagonal matrix.");	
					}
				}
				if (i == j) {
					d[i] = value;
				} else if (j - i == 1) {
					du[i] = value;
				} else if (i - j == 1) {
					dl[j] = value;
				}
				m_isFactorized = false;
			}
		};
		virtual property MatrixType Type {
			MatrixType get() { return MatrixType::DoubleGeneralTridiagonal; }
		};
		virtual property bool IsFactorized {
			bool get() { return m_isFactorized; }
			;}

		virtual void Invert() {
			throw gcnew LatoolException("Inversion of tridiagonal matrix is not supported by LAPACK.");
		};

		virtual void Factorize() {

			long n = m_rownum;
			long info;

			dgttrf_(&n, dl, d, du, du2, ipiv, &info);

			m_isFactorized = true;
		};

		virtual void Solve(IMatrix ^ b) {

			if (b->Type != MatrixType::DoubleGeneral) {
				throw gcnew LatoolException("Matrix b must be general matrix.");
			}
			DoubleGeneralMatrix ^ gb = (DoubleGeneralMatrix^) b;

			char trans = 'N';
			long n = m_rownum;
			long nrhs = gb->ColNum;
			long ldb = n;
			long info = 0;

			dgttrs_(&trans, &n, &nrhs, dl, d, du, du2, ipiv, gb->dat, &ldb, &info);
		};

		virtual void SolveWithSimpleDriver(IMatrix ^ b) {

			if (b->Type != MatrixType::DoubleGeneral) {
				throw gcnew LatoolException("Matrix b must be general matrix.");
			}

			DoubleGeneralMatrix ^ gb = (DoubleGeneralMatrix^) b;

			long n = m_rownum;
			long nrhs = gb->ColNum;
			
			long ldb = n;
			long info = 0;

			dgtsv_(&n, &nrhs, dl, d, du, gb->dat, &ldb, &info);

			m_isFactorized = true;

		};
	};
}